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Studying tlie jamming transition of granular and colloidal systems, has lead to a proliferation of 
theoretical and numerical results formulated in the language of the eigenspectrum of the dynamical 
matrix for these disordered system. Only recently however, these modes have been accessed experi- 
mentally in colloidal and granular media, computing the eigenmodes of the covariance matrix of the 
particle positions. At the same time new conceptual and methodological questions have appeared, 
regarding the interpretation of these results. In the present paper, we first give an overview of 
the theoretical framework which is appropriate to discuss the interpretation of the eigenmodes and 
eigenvalues of the correlation matrix in terms of the vibrational properties of these systems. We 
then illustrate several aspects of the statistical and data analysis techniques which are necessary 
to extract reliable results from experimental data. Concentrating on the case of hard sphere sim- 
ulations, colloidal and granular experiments, we discuss how to test the existence of a metastable 
state, the statistical independence of the sampling, the effect of the experimental resolution, and 
the harmonic hypothesis underlying the approach, highlighting both its promises and limitations. 
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I. INTRODUCTION 

Amorphous systems such as structural glasses, colloids, 
emulsions or granular matter still lack a satisfying de- 
scription. This is particularly apparent when one consid- 
ers the low temperature properties of glasses (T], Q or the 
sometime intriguing rheological properties of athermal 
amorphous systems, such as foams or granular mat- 
ter p, Q . Of particular importance is the understanding 
of what guarantees the mechanical stability of such sys- 
tems ■ The informations about the rigidity of a solid 
against collective particle motions are contained in the 
density of the vibrational modes D{uj): a system is sta- 
ble - at least linearly - if there are no unstable modes. In 
a continuous isotropic elastic medium, the invariance by 
translation implies that the vibrational modes are plane 
waves and the density of vibrational modes D{uj) follows 



the Debye law D{u 



'^0], where d is the space di- 



mension. By contrast, disordered solids exhibit common 
low-frequency vibrational properties that are completely 
unlike those of crystals. Disordered atomic or molecular 
solids generically exhibit a "boson peak", where n iany 
more modes appear than expected for sound |1CH13| |. 
Several empirical facts suggest that the presence of these 
excess modes is related to many of the original properties 
of amorphous solids [1, [13, [l^ . 

The jamming transition of frictionless granular mate- 
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rials lies at the threshold of mechanical stability, also 
known as the isostatic point. As a result of this equiva- 
lence, an excess of low- frequency modes vibrational arises 
and gives rise to a plateau in the density of states at low 
frequencies [l6| . This suggests that the zero-temperature 
jamming transition may provide a framework for under- 
standing at least some of the properties of amorphous 
solids. As a result of this an extensive theoretical and ex- 
perimental effort has been conducted towards the study 
of the vibrational properties of various amorphous sys- 
tems close to jamming. 

For frictionless repulsive soft spheres, the vibrational 
modes are obtained by diagonalizing the Hessian Kij = 

rP Y , defined as the matrix of the second derivatives 
of the pair interaction potential V{{ri\) with respect to 
the particles displacements around a given metastable 
state {r^} obtained numerically through conjugate gra- 
dient minimization. Above the jamming transition, it 
was shown that the onset frequency w* of the anomalous 
modes scales linearly with the distance from isostaticity, 
and that it is related to the existence of a diverging length 
scale in the response to an external perturbation|17l Il8l| . 
Beyond the frictionless, non-dissipative case, the den sity 
of states has also been studied for non-spherical |19l - [2]| 
and for frictionalp^. [23| soft particles. In the specific 
case of hard spheres below jamming, for which the po- 
tential is singular, an effective potential is first derived 
from the transfer of momentum during collisions in order 
to apply the same procedure [2^ . 

Experimentally, the density of states is traditionally 
accessed through light scattering techniques [2, 25, 26] or 
by computing the Fourier transform of the velocity auto- 
correlation function |27l | , as was recently done in thermal 



soft sphere simulations [28|] and in some granular experi- 
ments [23. Such methods do not give access to the spa- 
tial structure of the vibrational modes and only apply for 
thermally equilibrated systems. This is not the case of 
many systems of interest such as granular systems, foams 
and suspensions of large colloidal particles. For such sys- 
tems, recent advances in experimental techniques now 
allow tracking the real space dynamics of the individ- 
ual particles [30|. In principal one should thus be able 
to follow the same procedure as in simulations, namely 
average the particle positions in order to obtain a refer- 
ence state and then diagonalizing the Hessian around this 
state. However, the interaction potential V{{ri\) is usu- 
ally unknown. For granular systems or foams, the elastic 
part of the interaction is well understood but the pair- 
wise dissipation, either static friction or viscous damping, 
does not derive from a potential. In colloidal suspensions, 
the presence of hydrodynamics interaction prevent from 
deriving a simple two-body interaction rule. Also the ex- 
perimental resolution may be a limiting factor to perform 
such an analysis. 

An alternative approach is to follow the dynamics and 
perform a Principal Component Analysis (PCA) of the 
covariance matrix of the positions around the reference 
state of interest. This technique allows to extract the 



dominant components of the fluctuations - in the sense 
that they concentrate most of the correlated motion. It 
has been applied successfully in a range of fields, from the 
pathways of protein folding [sil, HH to financial mathe- 
matics j33| . Whatever the underlying dynamics are, the 
PCA acts as a filter which separates the correlated sig- 
nificant part of the fluctuations from the uncorrelatcd 
noise. 

For systems at thermal equilibrium, we shall recall that 
the eigenmodes of the covariance matrix of the positions 
are identical to those of the Hessian, and that the asso- 
ciated eigenvalues are related through a simple relation. 
Following this line several teams have recently attempted 
to extract the density of state from various experimental 
systems such as hard-sphere colloids [s^ - fSTj and vibrated 
granular materials [38|. While contributing to this ef- 
fort [13, [13, HI] , the authors of the present paper have 
realized that both the interpretation of the PCA analysis 
and its practical implementation requires some attention. 

In this paper, we would like to discuss both the promise 
and the limitations of this approach, as well as the set of 
hypotheses it depends on and the interpretation one can 
give to the experimental results. We pay particular at- 
tention to the influence of the real dynamics on the inter- 
pretation of the modes, and also point out the fundamen- 
tal differences between equilibrium and non-equilibrium 
systems. We discuss the differences between the idealized 
spring network to which the eigenfrequencies refer to and 
the actual physical system by developing the concept of 
shadow system^ first introduced in [3J] . We illustrate the 
discussion with simulations of jammed soft spheres obey- 
ing an over-damped Langevin dynamics. On the way, we 
will see that 

• (i) whenever the modes can be computed, they in- 
deed closely match those of the Hessian and are 
meaningful to the description of the dynamics, 

• (ii) the eigenspectrum of the covariance matrix of 
the position allows to compute the density of state 
only in the long time limit and if equipartition 
holds, 

• (iii) beyond the density of states, the covariance 
matrix of the positions is a rich source of informa- 
tion regarding the local rigidity properties of the 
system. 

We would also like to provide some practical guide- 
lines, which we hope will spare some time to those who 
would like to follow the above approach. We again illus- 
trate our purpose with several systems of interest. First, 
we use hard sphere simulations to address statistical is- 
sues, while experiments with granular systems close to 
jamming illustrate the possible anharmonicity of the dy- 
namics. Finally, two experiments with colloidal suspen- 
sions [H,!!^] illustrate convergence and resolution issues, 
in particular the problems that arise from an insufficient 
number of samplings [3^ . 
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At this point let us state that the present paper is 
written in a pedagogical spirit. It does not so much con- 
tain revolutionary ideas as instead the first reasonably 
complete and coherent treatment of the conceptual and 
practical issues surrounding this approach. We hope to 
have injected some wisdom on the use of PC A applied as 
to soft matter physics, and also of the broad relevance of 
the Marcenko-Pastur theorem. It is the result of many 
discussions with those colleagues of ours who pioneered 
this kind of delicate analysis, and we are grateful to them. 
It is our best wish that this paper provide an entry point 
to those entering the game. 



II. THEORETICAL BACKGROUND 

In this section we present rigorous derivations of the 
link between the covariance matrix of the positions and 
the Hessian matrix. Let's consider a system of N parti- 
cles following trajectories \r{t)) and let's assume one can 
define a statistical reference state |r°) = during 
a time interval T, where ( . ) denotes the average over 
time T. One then introduces the displacement vector 
\Sr{t)) = \r{t)) ~ |rO), the Hessian matrix K,j = 
defined as the matrix of the second derivatives of the pair 
interaction potential V(|r)) with respect to the particles 
displacements around the reference state jr"), and the 
covariance matrix of these displacements: 

Cp = (|Mt)) {6r{t)\). (1) 

Although the displacement vector \Sr{t)) depends on t, 
we will hereafter omit this dependence to simplify the 
notation. Cp is the empirical correlation matrix (i.e. on 
a given realization of the trajectories), that one must not 
confuse with the true correlation matrix of the under- 
lying statistical process. While in this first part we will 
assume that T is much larger than N ensuring the equiv- 
alence between time and ensemble averages and thereby 
the equality between the empirical and the true correla- 
tions, we shall see later that this is usually not the case 
in practical situations. Note that ensemble average here 
refers to the statistical average over the thermal fluctua- 
tions around a given glassy state. We will then have to 
refer to the Marcenko-Pastur theorem [s^ and its gen- 
eralizations to discuss the convergence properties of the 
spectrum of Cp to the one of the true correlation matrix. 

Cp is a real symmetric matrix, that can always be di- 
agonalizcd. We will see now how to relate its modes 
and eigenvalues to those of the Hessian, depending on 
the underlying microscopic dynamics. We first discuss 
two important classes of equilibrium dynamics, namely 
inertial Newtonian and fully over-damped Langevin dy- 
namics. The generalization to the case of an inertial and 
partially damped dynamics in the presence of a colored 
noise will allow us to discuss the case of athermal and 
mechanically excited systems. 



A. Equilibrium systems 

For systems in thermal equilibrium, statistical physics 
tells us that one can forget the dynamics and replace 
it by Gibbs statistics. One defines the partition func- 
tion Z (x J exp{— l3V)d\5r) and in the harmonic approx- 
imation V = ((5r|K|(5r)/2, so that it is straightforward 
to compute the correlation functions from the partition 
function. In particular, the correlation of the displace- 
ments reads 

Cp = kBTK-\ (2) 

Hence for equilibrium dynamics, the eigenvectors of 
Cp are simultaneously eigenvectors of K , while their 
eigenvalues are simply inversely proportional. We shall 
now illustrate how one can recover the above relationship 
in the cases of Newtonian dynamics and over-damped 
Langevin dynamics. The goal of this calculation is 
twofold. First it will point at the underlying hypoth- 
esis satisfied by equilibrium dynamics. Second it will 
provide a simple sketch of the more intricate calculation, 
which we will have to perform to discuss the case of non 
equilibrium mechanically excited systems. 



1. Newtonian dynamics 

In the harmonic approximation, that is linearizing the 
dynamics around the reference state jr^), one has for non 
dissipative dynamics 

\Sr) + '^\Sr) =0. (3) 

m 

This is the first and crucial approximation underlying 
the whole approach. Since the reference state is sup- 
posed to be mechanically stable, the eigenvalues Kg of 
Karc positive. Introducing D = K/m, also called the 
dynamical matrix, and where for simplicity we have 
supposed that all particles have a same mass m. The 
eigenmodes of D, also called the vibrational modes of 
the system, are then defined by D|Aq) = w^jAq), where 
LOq = -y/ Kq/m are the vibrational frequencies. 

The solution of Eq.Q is given by: 

\5r) =e~'^'\5r[0)), (4) 

where \5r{{))) is the displacement field at time zero. Re- 
placing this solution in Eq.([T|) and writing |(5r(0)) in the 
eigenbasis of D, {|Ag)}, one easily shows (see appendix 
IVII Al for derivation) that in the long time limit (long 
compared to the inverse of the minimal gap between ad- 
jacent frequencies), Cpis diagonal in the eigenbasis of 
D: 

Cp|A,) = al \\), (5) 
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where aq = {Xq\6r{0)) is the amphtudc of the initial 
condition on the mode \Xq). 

Now comes the third key ingredient: at equihbrium, the 
initial condition is thermalized and the energy is equally 
distributed among the modes. Each mode of frequency tUq 
and amplitude aq carries an energy maqUj'^/2 = kBT/2. 
Hence the eigenvalues Ag of Cp and the vibrational fre- 
quencies ujq of the dynamical matrix are related through 



0:„ 



(6) 



Practically, we have seen that for a Newtonian dynam- 
ics at equilibrium, the diagonalization of Cp does provide 
the vibrational modes of the system, the frequencies of 
which are proportional to the square root of the inverse 
of the eigenvalues of Cp . The largest eigenvalues of Cp , 
that is the most coherent motion corresponds to the soft- 
est modes of the system. 



that under such assumptions, Cp is diagonal in the basis 
of C, even for finite time. However it is only in the long 
time limit t » Tq for the largest Tq that one recovers 
the simpler expression 



A 



knT 



(10) 



where we have used the fluctuation-dissipation theorem 
r = 2^fcBT. 

Practically, we have seen that for an overdamped 
Langevin dynamics, the diagonalization of Cp again al- 
lows us to compute the vibrational modes of the system. 
The eigenvalues Xq are related to the relaxation time Tq 
of the dynamics in each mode, which are themselves re- 
lated to the oscillating frequencies Wg of the undamped 
system. 



B. Athermal systems 



2. Overdamped Langevin dynamics 

Again, in the harmonic approximation, one obtains the 
following Langevin equation for an over-damped dynam- 



\Sr{t)) = --\Sr) + -m), 



(7) 



where fi is the viscous damping and \ri{t)) is a white 
noise of amplitude F: {r]{t')\r]{t")) = TS{t' ~ t"), with 
F = 2fikBT since we consider equilibrium dynamics, for 
which the fluctuation-dissipation relation holds. 

To discuss this case, we introduce the operator £ = K//.J 
whose eigenvalue equation is C\Xq) = Kq/fi\Xq). 
The eigenvalues have the dimension of an inverse 
time, which one interprets as the relaxation time 



Tq = H/Kq 

mode I An 



ji/ (jii.ujT) of the system along the eigen- 



The solution of eq.([7]) can be written as 



\5r) 



-ct 



|5r(0)) 



M Jo 



-^(*-*')|77(t'))rft' (8) 



Following the same path as for the Newtonian dynamics, 
we show in the appendix IVH BI tliat Cp is diagonal in the 
eigenbasis of C : 



Cp\Xq) — 



2n' 



|A,), (9) 



where aq = {Xq\Sr{0)) is the amplitude of the initial con- 
dition on the mode |Ag). To derive the above relation, it 
is assumed that the initial conditions as well as the noise 
components on the modes arc both uncorrclated. Note 



For many systems of interest, such as grains, foams, 
suspensions of large particles, the thermal fluctuations 
arc orders of magnitude too weak to drive the dynamics. 
Such systems are thus generically out of equilibrium, the 
interactions lead to dissipation and some external forc- 
ing must be provided to drive the dynamics. The central 
thermal equilibrium relation equation ^ is based on as- 
suming both ergodicity of the fluctuations around the 
reference state during the observation period, and an un- 
derlying thermal canonical ensemble. Both assumptions 
are a priori violated in a non-equilibrium system. Also we 
have seen above that several hypothesis are necessary to 
derive a relation between the eigenvalues of Cp and those 
of the Hessian. Most of these hypotheses arc related to 
the correlation properties of the noise. These could be 
extended to non equilibrium situations. However, one 
key assumption is the equipartition of the energy on the 
modes. This property is very specific to equilibrium. It 
is well known that in out of equilibrium situations driven 
by a macroscopic forcing, the energy is in general not 
equi-distributed. On the contrary, the generic situation 
is that non linearities drive a cascade of energy from the 
large scale to the dissipative scale. For most systems, the 
spectrum of the energy fluctuations is not even known 
theoretically. 

Here we conduct the same kind of analysis as above 
for a homogeneously driven dissipative system. Vibrated 
grains experiments are typical realizations of this situ- 
ation. To generalize from the two equilibrium cases of 
Newtonian dynamics and Langevin dynamics, we now 
consider a non-equilibrium system with both inertia and 
damping and a colored noise spectrum. We choose to 
stick to a Langevin type description of the dynamics, 
where the damping is linear and single particle, as ap- 
propriate for particles in a newtonian fluid bath, but not 
necessarily for a large scale mechanical excitation. The 
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generic equations of motion for the linearized dynamics 
around a reference state |r*') is then: 

m\Sr{t)) + f^Mt)) = -^Mt)) + m) (11) 

where jj, is the viscous damping and |?7(i)) is a colored 
noise defined in the basis of the modes as 

with (r?g(t)%(t')> -r,V5(i-i')- (12) 

We assume that the noise correlations arc fully described 
by their first and second moments, i.e. a normal distri- 
bution. Here Tq is the amount of energy which flows 
into the system through non-equilibrium processes along 
this mode. It should be emphasized that in general this 
function is unknown and difficult to determine. 

Solving the dynamics projected on the modes of K , us- 
ing the same notation as for the Langevin case, K|Ag) ~ 
Hq\\q) (sec appendix I VII C|) . we find that Cpis diago- 
nal in the eigenbasis of K . At finite time one obtains 
complicated eigenvalues where the damping and the os- 
cillatory component of the dynamics interplay and one 
must distinguish the low ((/i/m)^ < AKq/m) versus the 
strong {{fi/m)^ > 4Kg/m) damping limits (equations l40l 
and HI] of the appendix). We have also assumed that 
noise and initial conditions do not cross-correlate, and 
that in ensemble-average, the initial conditions of differ- 
ent modes arc independent of each other. These condi- 
tions are in principle not so easy to satisfy if the noise 
is produced by a regular external excitation with a well- 
defined coupling to the modes. Fortunately, in the long 
time limit the system loses memory of its initial condi- 
tions and one recovers an expression similar to that of 
the equilibrium result 

CM^^W) (13) 

However one does not get rid of the violation of equipar- 
tition: the eigenvalues of Cp are the ratio of two ampli- 
tudes, one being the amount of energy the external forc- 
ing has put on the mode, the other being proportional to 
the mode stiffness. Hence the eigenvalues of Cp do not 
give access to the density of states of the system. 

III. INTERPRETATION 

We have seen in the above section that for systems in 
thermal equilibrium, one can in principle safely use the 
covariance matrix Cp of the positions around an equilib- 
rium state in order to obtain the stiffness matrix K , or 
equivalcntly the dynamical matrix D , and the eigcnfre- 
quencies of the system using the relations; 

K = keT Cp-^ (14) 



However, real systems can be vastly more complicated 
than the idealized Newtonian or Langevin equilibrium 
systems described above. Even a model system of 
jammed harmonic frictionlcss soft spheres obeying either 
Newtonian dynamics or treated with a conjugate gradient 
algorithm exhibits strong nonlinearities when approach- 
ing the jamming transition [iO, |4l|. The description of 
purely repulsive soft or hard spheres below jamming re- 
quires the introduction of an effective potential as defined 
by the transfer of momentum among the particles [l^l . In 
colloidal systems, electrostatic charges and the effect of 
solvent can significantly alter the pair interaction poten- 
tial including the development of long range interactions. 
Finally, we have seen that for out of equilibrium systems 
the possible violation of the equipartition relation pre- 
vents us from using relations ([H]) and ([15]). 

It is thus of primary importance to clarify the physical 
interpretation one can give to the eigenmodes and eigen- 
values of Cp . As we shall see below, there are two strate- 
gics. One may introduces a "shadow system" , which is 
by definition the thermally equilibrated system of which 
ksT Cp~^ is the rigidity matrix. As long as the hypothe- 
sis underlying relations (|14p and (|15p remains reasonable 
for the system of interest, this shadow system should in 
principle have the same eigenspectrum as the real one. 

If these hypotheses are not valid, in particular if the 
harmonic approximation or the equipartition are strongly 
violated, one should alternatively stick to Principal Com- 
ponent Analysis without expecting to relate the eigen- 
values of Cp to the vibrational properties. Note that the 
comparison with thermal systems of reference can still be 
done, but computing the spectrum of Cp from the den- 
sity of states of the thermal system instead of trying to 
compute the density of state of the athermal system from 
Cp . It can also be of interest, as we will illustrate be- 
low, to define, on the basis of the experimental data, a 
system of reference, the correlations of which have been 
suppressed. 

A. The "Shadow system": vibrational modes and 
density of states 

Let us define the shadow system explicitly. Suppose 
some experimental system of colloidal particles moving 
around a metastable configuration |r°). Once computed 
the effective spring constants fc°j^ between particles from 
the empirical inverse stiffness matrix K"*'' = kTC'^, 
the shadow system is defined as the Hamiltonian sys- 
tem of harmonic springs with stiffnesses k^^, strung be- 
tween particles of mass m sitting in the metastable con- 
figuration jr"). Note that the shadow system must not 
be confused with the model system the physicist has in 
mind when discussing his observations on the experimen- 
tal system. In the present case, despite possible compli- 
cated electrostatic and hydrodynamic interactions, the 
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FIG. 1: (color online) Left: Shadow system for a brownian system composed of soft harmonic spheres slightly above jamming 
at (f) — 0.86 and kT = 10~^. The contacts are plotted with widths proportional to their effective stiffness kcs, red with 
positive fcoff and blue with negative kcs- Center: Dependence of the effective stiffnesses k^s on the interparticle gap for different 
temperatures kT and — 0.86; compare to the step function for the harmonic potential set in the simulation. Right: Density 
of states D{u!) computed from the eigenspectrum of Cp for the same range of temperatures, compared to the one derived from 
the Hessian of the soft sphere potential. 



colloidal system is considered as a good experimental re- 
alization of over-damped soft spheres in a thermal bath 
- the model system the physicist has in mind. Whereas 
the model system obeys an ovcrdampcd dynamics, the 
shadow system obeys the Hamiltonian dynamics of per- 
fectly harmonic springs. 

For thermal systems at equilibrium, the spectrum of 
Cp can be interpreted as the vibrational modes of the 
shadow system and both the shadow system and the 
model system share the same density of states. How- 
ever one must remember that the vibrational properties 
of the model system are still different from those of the 
shadow system. First, even in the linear regime, it is well 
known that resonance frequencies of a system of springs 
are shifted by damping and eventually completely elimi- 
nated in the over-damped limit. Second, the model sys- 
tem may well escape the linear regime. This is what 
happens for frictionless disks close to jamming, because 
of the vanishing range of linear response, as was shown re- 
cently [4l| . In the following we illustrate this last remark 
with two computational "experimental systems" . The 
first system consists of poly-disperse packings of thermal 
frictionless soft spheres interacting through a harmonic 
potential, and simulated with over-damped Langevin dy- 
namics [4^ close to the jamming transition. The sec- 
ond one consists of thermal hard sphere packings simu- 
lated with event-driven Newtonian molecular dynamics 
just below the jamming transition. Initial states are ob- 
tained by shrinking simulated packings at jamming by a 
factor 6 € [10~^ — 10~^] which controls the distance from 
jamming [24| . 

Figure [T] displays the effective stiffnesses obtained for 
the soft sphere simulations slightly above jamming, at a 
very low temperature. We find effective stiffnesses consis- 
tent with the original connectivity and the original stiff- 
ness, which is a constant, here set to 1 for the harmonic 
potential. The only exceptions are occasional weak in- 



teractions with cage neighbors which are not contacts. 
As can be seen in Figure [T] center, at low temperatures, 
the stiffness as a function of overlap remains close to 
the T = step function. At low temperatures, as one 
would have expected for soft spheres above jamming, the 
shadow system is reasonably close to the "experimental" 
one. As a result the density of states D (w) computed 
from the Hessian matrix and the one computed from the 
spectrum of Cp closely match (see figure [3 right). 

When the temperature is increased, still remaining 
deeply within the glassy regime, we find that the stiff- 
ness as a function of the overlap broadens away from the 
T = step function, see Figure [1] bottom left. Equally, 
we find that the density of states obtained from Cp begins 
to differ substantially from the Hessian density of states. 
The larger temperatures allow the system to explore the 
potential beyond the range of linear response. 

The differences between the shadow system and the 
model system become strongly apparent if we cross the 
jamming transition at finite temperature. Figure [5] 
shows different aspects of the shadow system as a func- 
tion of packing fraction cj) at finite (low) temperature 
kT — 10~^. From the plot of the shadow system at 
(j> = 0.835, it is clear that the stiffnesses now deviate sig- 
nificantly from a step function. The gap distribution (top 
right) moves from distributions with only negative gaps 
(overlaps) above jamming, to a situation with only posi- 
tive gaps below jamming. Around jamming, we find both 
over- and underlaps within the same packing, and both 
contribute to the effective stiffness. Even though the 
stiffness-gap function has only slightly broadened from a 
step function, the different packing fractions sample dif- 
ferent parts of the curve (bottom left), and the resulting 
shadow system is significantly different from the model 
system. At high densities, we do recover a shadow system 
density of states which resembles the Hessian of the soft 
sphere potential (bottom right, compare to inset). At 
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FIG. 2: Top Left: Shadow system for a packing of thermal 
soft sphere just below jamming a,t (f> — 0.835 and kT = 10^®. 
The contacts are plotted with widths proportional to their 
effective stiffness fcefi, red with positive fceti and blue with 
negative kes- Top right: interparticle gap distribution for a 
range of packing fractions; with a negative mean gap above 
jamming and a positive gap below jamming. Bottom left: 
Dependence of the effective stiffnesses fcoff on the interparticle 
gap for different (f); the different densities sample different 
regions of the stiffness curve. Bottom right: Density of states 
D{ijj) computed from the eigenspectrum of Cp compared to 
the one derived from the Hessian of the soft sphere potential 
with an added log-potential (inset; see text). All curves for 
kT= 10"^ 



lower densities, the similarity breaks down, most clearly 
below jamming, where the soft sphere density of states 
has a large number of zero modes, and it is mechanically 
unstable, unlike the shadow system one. 

In the inset to the bottom right figure, we compare the 
shadow system density of states to the simplest effective 
model for thermal soft spheres: Veff = Vharm + kTViog, 
that is we have added an effective logarithmic potential 
for hard spheres to the soft sphere potential (see below 
for a disscussion of the log-potential). While Veff cap- 
tures part of the shift of the density of states to lower 
frequencies below jamming, it is not a good fit: all the 
shadow system density of states have a peak at interme- 
diate frequencies, unlike the density of states obtained 
through Veff- 

The effective log-potential we just used was first in- 
troduced for hard spheres slightly below jamming where 
the situation is again not so clear: the hard spheres po- 
tential is singular and, as for the soft spheres case, the 
system is a priori not mechanically stable anymore. On 
one hand one can derive an effective pair potential and 
effective stiffnesses by computing the forces amongst the 
particles via the exchange of momentum during the colli- 
sions [i^l • This potential scales like the logarithm of the 
interparticle gap. On the other hand, one can compute 
the effective stiffnesses derived from the inversion of Cr, . 




FIG. 3: Top: Shadow system for a packing of = 256 
thermal hard spheres just below jamming {(j) — 0.83878, 
= 0.83894). The contacts are plotted with widths pro- 
portional to their effective stiffness fccff, red with positive k^s 
and blue with negative feoff. Bottom Left: effective stiffnesses 
as a function of the gaps separating the particles: in black 
the effective stiffnesses derived from the effective logarithmic 
potential; in red those extracted from the inversion of Cp . 
Bottom right : Density of states D{ijj) computed from the 
eigenspectrum of Cp compared to the one derived from the 
Hessian of the effective potential. 



Figure [3] compares both effective stiffnesses for a system 
of = 256 hard disks just below jamming (</> = 0.83878, 
= 0.83894). We indeed find an effective connectivity 
of the packing even though in their mean positions, the 
particles are not touching. This is consistent with the 
effective potential approach and we also find k^s ~ ^/h^, 
where h is the interparticle gap, consistent with a loga- 
rithmic potential. Hence the shadow system remains a 
good approximation in this case too. As a consequence 
one can again compare the density of state D(uj) com- 
puted from Cp and the one derived from the Hessian of 
the effective potential. One observes on figure [3] that 
they match pretty well. Given the singular character of 
the hard sphere potential and the fact that the reference 
state is only a metastable state, this is a non trivial result 
which validates in a self consistent way both the approach 
in terms of an effective potential introduced in [24| and 
the use of Cp as a tool for extracting the density of states 
for slightly unjammed systems. 

Before closing this section, let us mention that we have 
followed the same analysis in the case of a real cxpcrimen- 



tal system composed of NIPA colloidal particles above 
jamming [isj . Conceptually, it is very similar to the soft 
sphere simulation we have just discussed. However, as 
already mentioned, the presence of a solvent and elec- 
trostatic interactions could alter the description in terms 
of a simple repulsive pair potential. We found that the 
effective stiffness is far above the noise for adjacent par- 
ticles only. This tells us that elastic interactions between 
neighboring particles dominate the data, as they should. 
Furthermore, the effective stiffnesses as a function of the 
interparticle separation could be measured and is com- 
patible with an harmonic or Herzian interaction. Hence 
also here the shadow system captures the essence of the 
vibrational properties of the real system under consider- 
ation. The reader interested by this experimental study 
should refer to the original paper. 



B. Principal component analysis 

So far we have discussed how to interpret the eigen- 
modes and eigenvalues of Cp for thermally equilibrated 
systems. We have seen in section Hi Bl that for non equi- 
librium systems, the distribution of the energy on each 
mode must be known before extracting the density of 
state. This in general is not the case. However the spec- 
trum of Cp and the associated modes still convey a lot of 
information about the dominant modes of relaxations. 

Such a path has been followed in analyzing different 
systems for which the underlying microscopic dynamics is 
unknown. In [ssj , the price fluctuations in financial mar- 
kets are analyzed by comparing the spectrum of eigen- 
values of the covariance matrix with the results of ran- 
dom matrix theory. The deviations from the purely ran- 
dom case contains the true information about the prices. 
In [sH the protein near-native motion is characterized. 
Although molecular dynamics simulations based on all- 
atom potential allows the identification of a protein func- 
tional motion on a wide range of timescales, the very 
large time scales are not easily accessible in simulations. 
A way to bypass this problem is to study the correla- 
tions of the fluctuations of the deviation of the protein 
backbone from its native configuration. The principal 
components of this matrix correspond to the collective 
modes of the protein which happen on larger time scales. 

Two of the present authors have used a similar ap- 
proach in the case of vibrated granular media (ssj . The 
experimental system consists in a 1:1 bidisperse mono- 
layer of brass cylinders. Mechanical energy is injected in 
the bulk of the system by horizontally vibrating the glass 
plate on which the grains stand and dissipated by solid 
friction. The experimental protocol produces very dense 
steady states with a packing fraction up to (j) — 0.8457. 
The stroboscopic motion of a set of 1500 grains is tracked 
in the center of the sample. The set-up, the quench proto- 
cols and the main properties of the system, are described 
in detail in [Ulil. 

We now briefly illustrate the methodology adapted 




10 



10" 



10" 







~ 1 

^ ~ ^ ^ -2 














•— • Experiment 
. . RM 

o 


-5/3 "■VZ/^'*"^^-*-,,^ 














10" 



10' 



10" 



FIG. 4: Eigenvalues of Cp for jammed grains under horizon- 
tal vibration above jamming. Top-left : Probability density 
of the rms displacement of each grain around it average posi- 
tion in the metastable state. Top-right : positive eigenvalues 
sorted in decreasing order. Bottom : probability density of 
the positive eigenvalues. On the center and right plots the 
black dots correspond to the real dynamics and the red square 
to the RMcr model; the dotted line correspond to the expected 
distribution for a crystal according to the Debye law 



in [38| to discuss the spectral properties of Cp . We have 
considered a subset of = 350 grains acquired during 
T = 100 time steps at a packing fraction $ — 0.844 
so that the metastable state of reference is well de- 
fined (no rattler, no structural rearrangement). Because 
r = dN/T, where d is the space dimension, is larger than 
one, strictly zero eigenvalues appear in the spectrum of 
Cp . We shall come back to this practical issue in the 
next section and restrict the discussion to the positive 
eigenvalues. The fluctuations of each grain i around its 
metastable position are Gaussian but the width, cr^, of 
the distribution is widely distributed amongst the parti- 
cles, according to the fat tail distribution p{a) ^ (j-i^+f^) 
plotted on flgure |4]-top-left , where /i w 6. 

The Marcenko-Pastur theorem, the purpose of which 
is to relate the spectrum computed at finite r to the one 
of the true correlation matrix could be used to obtain 
predictions but no explicit analytic form for the shape 
of the eigenvalue spectrum of Cp . Please see appendix 
IVIIDI for details. One alternative strategy is to define 
a Random Model of uncorrelated Gaussian variables the 
width (Tj of which are the experimental ones. We will 
refer to this model as the RMa case. The eigenmodes 
extracted from the real dynamics with eigenvalues larger 
than the largest one of the i?Mo- model contain the rele- 
vant information about the correlations. 

Figure |31-top-right, respectively flgure U-bottom, dis- 
plays the eigenvalues A™ sorted in decreasing order, re- 
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spectively the probability density p(A) extracted from 
the real dynamics as compared to those obtained for 
the RMa model. Note that there is a straightforward 
correspondence between these two curves since -^?ti(A) 
is exactly the inverse cumulated distribution of A and a 
power-law behavior A„i ~ in^^/°' translates into a density 
p(A) ^ A"'^"'^"'""^. The five largest eigenvalues for the ex- 
perimental system are clearly larger than the largest one 
obtained for the RM^ model. This comparison shows un- 
ambiguously that the top eigenvalues of Cp contain useful 
information about the dynamics of the system, and are 
not drowned in noise. It also demonstrates the existence 
of strong spatial correlations: by moving together, par- 
ticles achieve large collective fluctuations that would not 
develop otherwise. 

If one should not try to convert these eigenspectrum 
into density of states, it is however always possible to 
convert the expected density of states of a thermal crystal 
into the the spectrum of its dynamical matrix. Here, a 
2d crystal, with a density of states p{u)) ^ u'^^^ = u, 
with w cx A^/^, has p{X) ^ A~^, that is a = 1 and A,„ ~ 
m^^ as indicated on the figures by the green dotted line. 
The experimental data are consistent with an estimated 
a « 2/3, that is a slower decay of the spectrum than 
for the crystal case. One can interpret this result saying 
that there is a larger fraction of modes participating in 
the dynamics than in the crystalline case. 

One can also analyze the structure of the most sig- 
nificant modes. This was done in [ssj for this system 
of bi-disperse grains when approaching jamming from 
above, and it could be demonstrated quantitatively that 
there is a redistribution of spectral weight towards larger 
eigenvalues corresponding to softer modes and that the 
softest mode becomes more coherent and spatially orga- 
nized. However, as it was discussed above, for such a 
system where the energy is not distributed evenly over 
the modes, one cannot convert this information into the 
vibrational properties and the density of states of the 
system. The results presented here as well as in [s^ il- 
lustrate well how the PCA of the dynamics remains a 
fruitful tool of analysis without referring to the vibra- 
tional properties. 



IV. PRACTICAL GUIDELINE 

In the previous section, we have assumed that the data 
at hands were "perfect" in the sense that neither statisti- 
cal limitations nor resolution issues came into play. Also 
we did not consider the possibility of a strong anhar- 
monicity of the dynamics and we have not discussed the 
way, one should get rid of the rattlers, those particles 
which do not participate in the collective motion but in- 
stead "rattle" in their cages. 

Measuring Cp to a sufficient level of precision to ex- 
tract reliable information from it is not a simple affair. 
The goal of this section is to disentangle all the possible 
sources of practical problems that one faces when com- 



puting Cp and to propose methods and alternatives to 
obtain a faithful information from its spectrum. 

In subsection IIV A| we start by discussing a way to 
test the metastability of the reference state. Indeed all 
the approach is based on the fact the displacement field 
is purely composed of vibration around a steady refer- 
ence state. We illustrate this issue in the case of hard 
spheres, which by definition sit below jamming, where 
metastable states only have a finite lifetime. It is thus 
of crucial importance to control that the system does 
not escape the metastable state during the time window 
of the analysis. In subsection IIV Bl we then illustrate 
on the case of vibrated granular media, how one can ob- 
serve anhamonicity setting in when approaching the jam- 
mine transition from above. The reader may also refer 
to |4l| . where the same issue is discussed in the case of 
numerical simulations of soft spheres. As an immediate 
consequence of the above effects, the time window of the 
analysis may be seriously shortened as compared to the 
full duration of the numerical or real experiment. Hence 
statistical issues come into play. More specifically, we 
shall see in sections IIV CI and IIVDI how to deal with 
statistical independence of the data as well as conver- 
gence issues. In the case where the data sets have been 
obtained experimentally, resolution issues may also alter 
the analysis. We thus consider in subsection llV El a NIPA 
colloidal experiments to illustrate how convergence and 
experimental resolution interplay to affect the computa- 
tion of the density of states. Finally, identifying rattlers 
is usually a challenging matter. It emerges that calculat- 
ing the eigenspectrum itself is a good "filter" for rattlers: 
their motion appears only in a few low energy eigenval- 
ues which do not mix with the remainder of the motion. 
This is discussed in section IIVFI 

A. Metastability 

Before starting any of the analysis discussed in the pre- 
vious section, one must have checked that the dynamics 
is purely composed of fiuctuations around a well defined 
reference state. This can be done both a priori, in order 
to select a good time window to perform the analysis, 
and a posteriori on the basis of the properties of Cp . 
A natural tool to characterize the dynamics is the self- 
intermediate correlation function, defined as 

C,{to,t) = ( cos{q .Ar,{t,to) )j, (16) 

where the average is made on the particles but not on the 
initial time. Ari(t,to) = ^j{t)~^j{to) is the displacement 
of the particle j between to and t + to, q is a vector 
whose amplitude is given by q = tt/o and a is a length 
scale which sets the amplitude of the displacement above 
which the particle motion induces a change of metastable 
state. There is obviously a large part of arbitrariness in 
such a definition and traditionally a is set to a fraction of 
the particles diameter, so that a metastable state more 
or less corresponds to a given configuration of neighbors. 
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FIG. 5: Probing metastability. Top : Cq{to,t) for a system 
of hard spheres below jamming (TV = 1024 and (j> = 0.841947; 
(j>j = 0.841959). The time unit is tf, the time interval at 
which particle positions are stored, and corresponds to 800 
collisions per particle. Each interval delimited by dashed lines 
has 200 snapshots, during which we compute Cp and its basis 
Bi of eigenvectors {|Am)}. Bottom : Robustness of the first 20 
modes as defined in Eq. (|17|) between pairs of intervals Bi x Bj 
as indicated in the legend. The dashed line indicates the value 
of R(m) for purely random eigenvectors. 

Figure 0-top displays Cq(to,t) computed for a system 
of iV = 1024 bidisperse hard discs the dynamics of which 
has been simulated using molecular dynamics sli ghtl y be- 
low jamming {(j) = 0.841947; cj).; = 0.841959) [li, El. 
Here a has been set to the radius of the smaller particles. 
One observes very well defined plateaus, during which 
the system is trapped in a metastable state, separated 
by a quick relaxation event that we call crack. Typi- 
cally, the size of the plateau defines T, the interval of 
time available to compute Cp . There are however some 
cases where one cannot identify the plateaus so easily. As 
a matter of fact this will happen each time the system 
size is large as compared to the typical size of a crack. 
When this is the case, one can always cut the system into 
smaller subsystems and do the analysis in each subsystem 
as illustrated in [s^ ■ 

As long as the system remains in the same state of 
reference, the basis of eigenvectors should essentially re- 
main unchanged. This can be checked by dividing the 
time interval of duration T in smaller subintervals and 
computing an indicator of the modes robustness, defined 
as [38j : 

j=+M 

Rim) = (17) 

j = ~M 

where {^m\^'m+j) is the scalar product between the 
modes computed during two successive observation win- 
dows of duration T' < T. Setting M to a small but 



strictly positive value allows neighboring modes to possi- 
bly exchange their rank. As long as the largest eigen- 
values are concerned, we could check that the results 
are basically independent of M for M > 2 and we fixed 
M = 2. Figure [5]-bottom displays R{m) for the twenty 
largest eigenvalues computed between intervals defined 
as indicated on figure [S]-top by vertical dashed lines. The 
robustness of the the first modes is higher when com- 
puted amongst intervals belonging to the same plateau. 
As soon as one compute the modes on an interval in- 
cluding the crack, the robustness decreases sharply and 
rapidly reaches the baseline level R{m) = {2M + 1)/2N, 
one would obtained for purely random vectors. Note also 
that the robustness is not restored when the two intervals 
belong to different plateaus indicating that the system 
has indeed relaxed from one metastable state to another. 



B. Harmonicity of the dynamics 

In the previous section, we have identified a metastable 
state and checked a posteriori that the basis of eigen- 
vectors remains unchanged during the lifetime of this 
state. As a matter of fact, one must also check that 
the dynamics around this state is purely harmonic. For 
an equilibrium system close to jamming, the structure 
of the packing is frozen and fluctuations of the parti- 
cles are thermal, so that assuming Gaussian fluctuations 
around the metastable state sounds reasonable. However 
it was recently claimed that repulsive contact interac- 
tions make jammed particles systems inherently anhar- 
monic (4l| . The authors argue that at very low temper- 
ature the breaking and forming of inter-particles contact 
is responsible for anharmonicity in the sense that the re- 
sponse does not remain confine to the original mode of 
excitation. Such a definition of harmonicity is very strict, 
and how the result depends on system size, temperature 
and distance to jamming is still a matter of debate. It is 
clear that for mechanically excited athermal and dissipa- 
tive systems, such as shaken macroscopic grains, there is 
no reason why the dynamics should be harmonic. 

While studying horizontally shaken grains, two of the 
present authors proposed to check the harmonicity of the 
dynamics in a naive but experimentally accessible way, 
namely by computing the distribution of the individual 
fluctuations: p{5xi/ai), where 5xi is the position fluc- 
tuation along a given direction and Ui is the root mean 
square displacement of particle i. Figure [5] presents such 
distributions obtained for a system of iV = 1550 brass 
disks shaken horizontally on a oscillating plate just above 
jamming = 0.8417, see [l^ for details on the ex- 
perimental set up and protocol). The distributions are 
computed for four different packing fractions cj) and four 
durations T of the window of observation. They are then 
ensemble averaged over the lO^/T intervals provided by 
the full dataset. 

The distributions shown in figure [5] highlight some im- 
portant characteristics of the dynamics. The parameter 
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FIG. 6: Distributions of the position fluctuations p{Sxi/ai) 
for four values of the packing fraction and four durations of 
observation T. We focus on the top of the distribution, which 
is compared to a Gaussian (continuous black line) . The green 
hatched line separates the distributions well described by a 
Gaussian from those deviating from this harmonic description 



space {(/), T) can be divided into two regions, as illustrated 
by the green hatched line: for small enough observation 
duration T or large enough packing fractions, the dis- 
tributions are unimodal with a Gaussian core: particles 
jiggle around a well defined average position; for longer 
T or smaller densities, the distribution starts developing 
a flat top, with a poorly defined maximum. This sug- 
gest that on these longer observation times, the average 
position of a significant part of the particles is not well de- 
fined anymore. Particles either drift slowly or even find 
(collectively) another metastable position, as suggested 
by the double peak observed in the case (j) = 0.8417 
and T = 10^, i.e. for the loosest packing fraction and 
the longest observation time. This means that over long 
time scales, the evolution of the average position becomes 
comparable or even larger than the fluctuations, and it 
becomes meaningless to describe the system in terms of 
small vibrations around a fixed metastable state. For an 
infinite size system, some rearrangement always happens 
somewhere, and the covariance matrix Cp is always ill- 
defined. The "allowed" time scale Tniax(0, L) is expected 
to scale inversely with the system size L. 



C. Statistical independence 

Once a reference state is identified during a time win- 
dow of duration T, one needs to know whether the parti- 
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FIG. 7: Testing statistical independence during a given 
metastable state for a system oi N — 256 hard disks below 
jamming {(j> = 0.83878; = 0.838865). Top : sketch of the 
time discretization procedure (see text for details). Bottom 
: {R{m)) for the three discretization scheme as indicated by 
the legend. 



cle positions stored on successive snapshots are indepen- 
dent, or at least sufficiently uncorrelated to compute Cp . 
Again computing the robustness of the basis of eigen- 
modes is a good way to check the validity of the analysis: 
if the particle positions are not independent the correla- 
tor Cp is badly estimated and the eigenmodes are mostly 
composed of noise. To illustrate this point we use again 
the hard discs system just below jamming (0 = 0.83878; 
(j>j = 0.838865), now with N = 256, for which we could 
identify a long lasting metastable state the total dura- 
tion of which, T, corresponds to 1.6 x 10^ collisions per 
particles. 

Figure [7]-top sketches the procedure. The total acqui- 
sition is divided in 10 intervals. Each of them contains 
800 snapshots of the positions, separated by t/ = 200 col- 
lisions per particle. Cp is computed and diagonalized on 
each interval, and the average robustness among succes- 
sive intervals (i?(m)) is evaluated for the first 20 largest 
eigenvalues. To preserve the same statistics, together 
with a reduced interval of time between the snapshots, 
we then subdivide each interval in 10 subintervals, each of 
them containing again 800 snapshots of the positions, but 
now separated hy tf = 20 collisions per particle. And the 
procedure is iterated once more leading to 800 snapshots 
of the positions separated hy tf =2 collisions per parti- 
cle in each interval. Figure [7]-bottom compares {R{m)) 
for the three successive level of time discretization. Al- 
though the system remains in the same metastable state, 
the number of snapshots and the averaging used to com- 
pute Cp are identical in all three cases, one clearly ob- 
serves the convergence of (i?(m)) towards larger values 
when the snapshots are well separated. In the present 
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case the finest discretization is obviously irrelevant and 
a separation of successive snapshots larger then about 50 
collisions per particle is necessary. 

Altogether, whether it comes from experimental lim- 
itations of the acquisition rate or from a microscopic 
timescale inherent to the dynamics, there is always a 
minimal timescale separating the accessible and useful 
snapshots of the particle positions. Since we have also 
seen that the analysis is confined to the lifetime of the 
reference state, we are now in the position to face the 
convergence issue of the cigcnspectrum as discussed by 
the Marccnko Pastur theorem. 



D. Convergence of the spectrum (small sample 
limit) 

The Marcenko-Pastur (M-P) theorem [H, [H IH, 
applies to the ensemble of Wishart random matrices; 
that is the ensemble of p x p-matrices constructed by 
the addition of n independent samples XiXj, where the 
Xi (i € l-.-p) are distributed according to either a multi- 
variate Gaussian distribution or a distribution with finite 
higher moments, similar to the conditions of the central 
limit theorem. The M-P theorem then predicts the eigen- 
value distribution of such matrices in the limit n — > oo, 
p — >■ oo, r = p/n finite. However, as one can see in 
appendix IVII Dl there is no explicit form for the distri- 
bution. Also, the theorem tells nothing about the con- 
vergence of the eigenmodes. This is why we concentrate 
here on the numerical characterization of the convergence 
in a specific and well controlled situation. This will al- 
low us to discuss both the shape of the spectrum and 
the relevance of the modes themselves. Also, it will allow 
us to examine the case where r > 1, when a finite part 
of the spectrum is strictly zero. The reader interested 
in analytical predictions may refer to the original papers 
[H m m, 113: as well as to appendix EUdI 

We consider the same system as above, with N ~ 256 
hard disks just below jamming (0 = 0.83878; ~ 
0.838865) within a long lasting metastable state of to- 
tal duration T = 4.6 10^ collisions per particles and we 
choose to retain as independent particle positions sepa- 
rated by 200 colhsions per particle, so that n = 23000. 
In section IIII A[ we have shown that the spectrum of 
Cp and that of the effective dynamical matrix match 
when the whole set of data is used to compute Cp . Here 
we discuss how much the spectrum and the modes of 
Cp deviate from their asymptotic behavior when the time 
interval used to compute Cp is artificially reduced such 
that r = 2N/n varies in the range [0.021 - 2.28]. 



Figure [5] displays the spectral properties of Cp for in- 
creasing values of r. As soon as r > 0.5, the smallest 
eigenvalues are seriously underestimated, respectively the 
largest frequencies are overestimated. For r > 1, strictly 
zero eigenvalues replace the lowest eigenvalues. On the 
contrary, the largest eigenvalues, respectively the low- 
est cigenfrequencies, arc surprisingly robust, even in the 
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FIG. 8: Convergence of the spectrum of eigenvalues of Cp for 
a system of A'' = 256 hard disks below jamming (0 — 0.83878; 
(f)j — 0.838865) for decreasing values of r = 2N/n as indi- 
cated in the legend. Top: Eigenvalues (left) and correspond- 
ing eigenfrequencies (right) sorted in decreasing, respectively 
increasing order. Bottom: Eigenvalues spectra (left) and cor- 
responding density of states (right). 
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FIG. 9: Matrices of the scalar product between the modes 
of Cp and the modes of the dynamical matrix for a system 
of A = 256 hard disks below jamming {(p = 0.83878; (pj = 
0.838865). The rows are the indices of the modes in increasing 
order of their frequency and the colorbar indicates the value 
of the scalar product. Left: r = 2N/n = 0.021 Right: r = 
2N/n = 2.78. 



cases r > 1. Because of this robustness, when one is 
interested in the low frequency part of the density of 
states, it is recommended to normalize the density such 
that the area under the curve equals the fraction of non- 
zero eigenvalues. As a result, the density of states D{oj) 
conserves its shape as long as r < 1.5 and only for larger 
r the large frequencies part of the density of states starts 
to be significantly depleted. 

Finally figure H] illustrates the robustness of the modes 
themselves. To do so, we define the matrix of scalar prod- 
uct between the two eigenbasis: the modes of Cp and the 
modes of the effective dynamical matrix. If the modes 
were strictly the same, this matrix should be the iden- 
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FIG. 10: Density of states for a system of hard disks be- 
low jamming (0 = 0.841947; 4)j = 0.841959) extracted from 
Cp either computed using the whole system of A'' = 1024 hard 
disks (in blue) or cutting the system in 4 sub-systems with 
j^sub _ each (in red) as compared to the density of states 
obtained from the effective dynamical matrix computed for 
the whole system (in black). Inset: zoom on the lowest fre- 
quency part of the density of states. 



tity matrix. Figure [9]-left shows the matrix of the scalar 
products for r = 0.021, from where we observe that the 
modes of Cp project on a very small number of modes of 
the effective dynamical matrix, as already pointed at in 
section ITlI Al We see that even for r = 2.78, the modes 
of Cp associated with the largest eigenvalues still project 
on a small number of modes of the effective dynamical 
matrix, while the modes of Cp with a higher frequency 
are spread across the whole spectrum of the effective dy- 
namical matrix. 

Altogether, as long as r < 1.5, the most significant 
eigenmodes and eigenvalues of the dynamic matrix as 
well as the qualitative features of the density of states are 
well captured by the analysis of Cp . Practically speak- 
ing, when very long lasting metastable states are avail- 
able it is thus often advantageous to average the spectrum 
or the density of states on successive time windows, pro- 
vided that the duration of each of them satisfies r sa 1 . 

When there is no enough snapshots to recover the 
whole spectrum of modes, an alternative strategy is to 
reduce the system size in order to increase r. In prin- 
ciple one can divide the system in smaller subsystems, 
and compute the density of states in each of them. This 
strategy allows to better access the low frequency part of 
the density of states, although the lowest accessible fre- 
quency increases as the inverse of the system size. Fig- 
ure [To] and figure [11] illustrate the application of such a 
strategy to a system of iV = 1024 particles below jam- 
ming (0 = 0.841947; (pj = 0.841959). In this example, 
the metastable state has a total duration of T = 1.6 x 10^ 
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FIG. 11: Matrices of the scalar products between the modes 
of Cp and the modes of the dynamical matrix for a sys- 
tem of N — 1024 hard disks below jamming. Left: Modes 
of Cp computed on the whole system. Right: Modes of 
Cp reconstructed from the modes computed separately in four 
sub-systems of size N^ub — 256. The horizontal (resp. the 
vertical) indices label the modes of the dynamical (reps, the 
covariance) matrix. The color bar indicates the value of the 
scalar product. 



collisions per particles and we store independent particle 
positions separated by 80 collisions per particle. This 
means that n = 2000 snapshots and r = 1.02. From 
figure H] we know that for this value of r the large fre- 
quencies are underestimated. By cutting the system in 
4 boxes, each one has r = 0.25 and we are now able to 
recover the whole spectrum. In figure [10] we compare the 
density of states D{u;) extracted from Cp either using the 
whole system of = 1024 hard disks (r = 1.02) or cut- 
ting the system in 4 sub-systems with N^^^ — 256 each 
(r = 0.25) to the density of states obtained from the ef- 
fective dynamical matrix computed for the whole system. 
Whereas the spectrum obtained from Cp computed with 
the whole system only captures the very low frequencies, 
the one computed for the sub-systems matches the spec- 
trum from the dynamical matrix on the whole range of 
frequencies. 

Concerning the robustness of the modes, once the 
modes have been computed in each subsystem, we first 
reconstruct artificial modes for the whole system by past- 
ing together the modes of each subsystem. Note that the 
number of modes in these 2 basis is not the same: for 
the whole system there are 2A^p = 2048 modes while 
in each subsystem has 2N^'^^ = 512 modes. Figure [11] 
shows the matrices of the scalar products between the 
modes of the dynamical matrix and those of Cp , either 
computed on the whole system (left) or reconstructed as 
described above from the modes computed on the subsys- 
tems (right). Despite the artificial nature of the modes 
obtained from the computation of Cp in each subsystem, 
their projection on the modes of the dynamical matrix 
remains rather concentrated, especially when one remem- 
bers that each mode of Cp in this case projects approxi- 
mately on four modes of the dynamical matrix. 

Let us end this section by quoting a recent work [i^ 
which shows that in general the modes are modified by 
the truncation of the full real-space correlations and pro- 
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poses windowing strategies to best avoid these boundary 
conditions artifacts. Implementing such strategies on top 
of the segmentation proposed here would certainly de- 
serve proper attention. 

E. Convergence and experimental resolution (large 
sample limit) 

Even in the most favorable case of all experimental sit- 
uations discussed above, we have not yet addressed the 
issue of the experimental resolution. As detailed in ap- 
pendix |VITd1 the Marcenko-Pastur theorem allows us to 
compute the impact of a finite resolution on the spectral 
properties of Cp at a finite r = 2N/n, where n is again 
the number of independent snapshots of the particle po- 
sitions. 

To do so we assume that the applicability conditions 
of the theorem detailed in appendix IVII Dl are met. The 
finite moments condition corresponds broadly speaking 
to a long-lived metastable state with short-time fluctua- 
tions of the particles in well-defined cages. The metasta- 
bility condition can be tested through the methods of the 
section IIV Al For an equilibrium system, the short-time 
fluctuations of the particles in their cages are thermal and 
can hence well be approximated as gaussian within the 
range of linear response. For a non-equilibrium system 
such as shaken granular systems, this condition needs to 
be verified experimentally. The independence of the dif- 
ferent samplings of the system can be tested through the 
methods of section ITV CI 

If we have now constructed Cp in these conditions, how 
does the experimental resolution affect the properties of 
the matrix and how does it interplay with the finitencss of 
r? How docs it modify the spectral properties of the true 
correlation matrix, the one would obtain in the limits of 
infinite resolution and r — >■ 0? We model the resolution 
as an additive Gaussian noise of variance e^, where e is 
the resolution-induced uncertainty of the particle posi- 
tions. If Co is the true correlation matrix, i.e. the one of 
the shadow system, one easily shows that for Gaussian 
fluctuations of the particle positions, the finite resolu- 
tion correlation matrix is Cq = Cq + e^I, where / is the 
identity matrix. Then the measured Cp is Cr, such that 
lim„_j.ooC'r = Cq. The Marccnko Pastur theorem is still 
valid and one finally relates the eigenvalues distribution 
of Cr and the true correlation matrix Co- 

As for the case without noise, there is no explicit form 
for the distributions, but we show in appendix IVII D I that 
the mean and variance of the eigenvalue distribution of 
Cp converge with r — > and e//io — > in the following 
way: 

yUr = Mo + 

= al + r^il + re^ (18) 

where the quantities labeled refer to the exact eigen- 
value distribution and the ones labeled r to the mea- 
sured quantities. Note that for uncorrclated particles /ig. 



would be the typical cage size squared, i.e. the value 
of the plateau in the mean square displacement curves, 
preceding the onset of diffusive behavior. Interestingly, 
while the resolution directly impacts the mean, it only 
changes the width at ©(e**), and the correction is more- 
over multiplied by r. As a result the noise dependence of 
the eigenvalue distribution is in fact much weaker than 
one might have expected. 

We now present the results of a statistical model specif- 
ically developed to test the validity of the experimental 
results for colloidal NIPA particles [l^. In the experi- 
ment. N = 3600 particles forming quasi two dimensional 
packings deep in the glassy phase were imaged with a con- 
focal microscope, yielding n = 30000 independent sam- 
ples of their positions per run. In the language of the 
Marcenko-Pastur theorem, this corresponds to r = 0.24. 
The combination of optical resolution and sub-pixel ac- 
curacy particle tracking led to an estimated uncertainty 
in the particle positions of e = O.OOT^m. 

We build a 'random vector' model by constructing a 
model Cp where we repeatedly sample from an uncor- 
rclated multivariate normal distribution coresponding to 
the experimentally measured cages sizes, and with added 
noise at the experimental amplitude. The covariance ma- 
trix of the model is then Cij = {\i + t^)5ij, where the 
are the diagonal elements of one of the measured Cp from 
. While this model resembles the 'Random Model' in- 
troduced in section IIII Bl detailing PCA, our aim here is 
different: we only intend to illustrate equation (TS] with a 
model where we know the exact eigenvalue distribution 
in the absence of noise and for r — > 0. By construction, 
this is just P{Xi), the distribution of diagonal elements 
or cage sizes which enter the model. 

Figure [12] (top-left) shows the convergence of the nu- 
merical eigenvalue distribution Pr{X) as r — > at fixed 
noise. The eigenvalue distribution corresponding to the 
experimental value r ~ 0.24 is shown in black. It appears 
to be relatively close to the exact distribution shown in 
pink. The numerically obtained first two moments of the 
distributions exactly match the theoretical result from 
equation (fT8|) (fig. [T2Ktop- right). The curves ndicate that 
for the experimental value of r (green vertical line), the 
error made compared to the true converged noiseless dis- 
tribution is of the order of 10%. 

For equilibrium systems, one is mostly interested in the 
density of states. Unfortunately, when performing the 
change of variables w = ^/ksT /raX, there is no straight- 
forward and reliable way to calculate the moments of 
the distribution and we have to rely mostly on numer- 
ical results. Only the asymptotic values in the limit 
r ^ could be obtained (see appendix IVII Dl) . In the 
bottom row of figure [12] wc have applied the transfor- 
mation uj — 1/ -s/A numerically to obtain the density of 
states and its first two moments for the 'random vector 
model. The density of states obtained for the experi- 
mental value of r = 0.24 again matches pretty well the 
asymptotic one. Note that we find a fortitous compen- 
sation in the moments of the distributions: in the limit 
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FIG. 12: Statistical convergence analysis for the NIPA system [l^, using the 'random vector' model (see text). Top left: 
Numerical eigenvalue distribution for different values of r as indicated in the legend. The distribution associated with the 
experimental value r = 0.24 is plotted in black and the asymptotic r — > one is plotted in pink. Bottom left: Same for the 
density of states. Top right: Numerical (dots, subscript 'sim') and theoretical (dashed lines, subscript 'th') results for the mean 
and variance of P(A) normalized by the exact result. The experimental value of r = 0.24 is indicated by a green line. Bottom 
right: Same for the density of states, numerical result only (dots and dashed lines, subscript 'sim'). The green line marked '30 
000 frames' corresponds to the number of samples in the experimental NIPA system [l3| . 



r — 0, one can show that a finite resolution systemati- 
cally leads to an underestimation of both the mean and 
the variance. Since working at finite r always leads to 
an overestimation of these moments, it compensates the 
above underestimation. In the present case, the experi- 
mental value of r (green line) is very close to the point 
where the compensation is perfect. Note that this is pure 
coincidence related to the specific values of the moments 
of the true distribution, and that even so it does not val- 
idate the detailed shape of the distribution. 

In the case of the experimental data obtained with 
NIPA particles [l^l, and for which correlations arc 
present, the authors tested convergence of the density 
of states based on the methods outlined here and found 
scaling of mean and width consistent with our example. 
The error from experimental and statistical resolution 
was estimated in the 10% range at r = 0.24. 

We can make a similar comparison for the hard sphere 
colloid setup of Ghosh et al. [3^. For this experiment, the 
estimated number of particles is = 2000, and approx- 
imately n = 4000 samples at a resolution of e = 0.03/im 
were taken to obtain the DOS. While the resolution ef- 
fect is comparable to the NIPA setup, we obtain r « 1. 



This introduces a large error into the measurement of the 
DOS, similar to the r = 0.95 curves in the left panels of 
Figure [m We estimate an error of 40% and 200% on the 
mean and the width of the DOS, respectively. However, 
as detailed in the previous section, even at this r it is still 
possible to obtain the relevant low energy modes [49| . 



F. Rattlers 

Close to jamming, there are few particles, the so called 
rattlers, which do no participate to the rigidity of the 
structure: removing them does not alter mechanical sta- 
bility. When computing for instance the average coordi- 
nation of the packing, these rattlers need to be identified 
and excluded. For simulations of soft spheres, the proce- 
dure is straightforward. The number of contact is known 
for each particle and those which have less than two con- 
tacts are considered as rattlers. For simulations of hard 
spheres just below jamming, the contacts are defined by 
computing the transfer of momentum during successive 
collisions. Again the particles with less than two con- 
tacts arc defined as rattlers, but also those which have 
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FIG. 13: Eigenvalues of Cp and participation ratio as defined 
in the text for a system of A'' = 256 hard disks below jamming 
{4> = 0.83878; 0j = 0.838865) with its rattlers (black curves) 
and having removed them according to difi'erent values of the 
thresholding criteria ^ as indicated in legend. See text for the 
definition of ^. 



a significant smaller momentum transfer per contact as 
compared to the average. Such a criterion was proposed 
in [50[, where it was shown that particles with momen- 
tum transfer per contact smaller than 2% of the average 
value could safely be considered as rattlers. Varying this 
threshold in the range 0.5 — 5% did not alter the results. 

We now offer an alternative way to identify rattlers in 
experiments by using the Principal Component Analysis 
(PCA) itself and compare it to the above method used 
for hard spheres. As a matter of fact, the PCA diago- 
nalizes the dynamics on its dominant features. Since the 
rattlers are not tightly confined by their neighbors, they 
move significantly more than other particles. Also such 
motion are by definition localized on a single particle. 
Hence rattlers correspond to very singular kind of exci- 
tation and thereby tend to jeopardize the modes struc- 
ture. Diagonalizing Cp and analyzing the modes with 
high eigenvalues and a very high participation ratio is 
thus a good strategy to identify rattlers. Here the par- 
ticipation ratio of a given mode m is simply defined as 
Pirn-) = X]fc=i('^''D^' where (5r™ is the displacement of 
the particle i on the mode m. Given that the modes are 
normalized, '^f^iiSrf^)'^ = 1 so that a perfectly dclo- 
calized mode has a participation ratio of 1/A^, whereas 
a mode fully localized on a single particle would have a 
participation ratio of order one. 

The method is illustrated in the case of a system of 
N = 256 hard disks just below jamming (0 = 0.83878; 

= 0.838865) within a long lasting metastable state of 
total duration T = 4.6 x 10^ collisions per particles with 
n = 23000 independent snapshots and r ~ 0.021. Fig- 
ure [13] shows the spectrum and the participation ratio 
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FIG. 14: Eigen modes of Cp for a system oi N = 256 hard 
disks below jamming (0 = 0.83878; (jij = 0.838865) with its 
rattlers (top row) and having removed them (bottom row). 



of the modes of Cp computed with and without rattlers. 
One easily identifies the modes with a simultaneously 
large eigenvalue and high participation ratio. Three of 
these modes are plotted on figure [T4l-top. Here also, one 
immediately identifies the rattlers with huge displace- 
ments as compared to the other particles. The rattlers 
in each modes are then defined as the particles having a 
displacement amplitude larger than ^ times the average 
displacement inside the mode: i5r™ > {Sr"^) +^ag^, with 
cr^ being the standard deviations of the displacement on 
the mode m. One must fix a threshold; however because 
the modes selected by the procedure concentrate most of 
the motion on the rattlers, it is very easy to fix the thresh- 
old and the results are robust to the choice of as it can 
verified in the figure[T3l Once the rattlers are identified - 
sometime in several modes one eliminates them and re- 
compute Cp . One then easily checks that the extremely 
localized modes have disappeared, together with the rat- 
tlers, figure [T4l-bottom. Note that most of the particles 
identified as rattlers using this criterion are identical to 
those identified with the previous criterion used in [soj . 
The mismatch between both criteria is of the order of 
10% of the rattlers identified. 



V. CONCLUSIONS 

We have reviewed various situations for which the 
study of Cp , the covariance matrix of the positions of 
a system of particles can be used to infer dynamical 
properties of the system. It has been emphasized that 
for systems without equipartition of the energy, or with- 
out an independent knowledge of the energy repartition 
amongst the modes, the spectrum of Cp cannot be trans- 
posed into the density of states of the system. Still, it 
conveys relevant informations about the dynamics and 
physicists would gain in exploring the path already fol- 
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lowed. We have also insist on the physical interpretation 
of the modes obtained via the analysis in terms of the 
vibrational properties of the shadow system, which must 
not be confused with the experimental system, even an 
idealized one. 

We then have seen the numerous steps one has to go 
through in implementing the analysis. Metastability of 
the reference state, harmonicity of the dynamics, statis- 
tical independence, convergence and experimental reso- 
lutions are all key elements to care of. How to deal with 
rattlers was finally discussed. Following the whole pro- 
cedure requires as always a bit of physical insights, a 
good knowledge of the system of interest and patience. 
We have not discussed the new physics one can finally 
extract from the analysis. It certainly depends on the 
specificity of each system, but in any case one should 
remember that the analysis remains linear in essence. 

We hope that in this light the present review of the 
methods associated to the correlation matrix approach 
will be useful to the community. 
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VII. APPENDIX 



where \Sr{0)) is the displacement field at time zero. Re- 
placing this solution in Eq.([T|), one has 

Cp = ( e^^VKO))(^r(0)|e*^*) ; (21) 

then writing \dr{0)) in the eigenbasis of D , {|Ag)}: 

|MO))=E"'?l^'?)> (22) 
q 

where aq = {Xq\dr{0)) , is the amplitude of the initial 



condition on the mode \Xq), one obtains 




a,a^e-*(-'-'')*|A,)(Afc| 



(23) 



Provided that the time average (.) is performed on a large 
enough interval, large as compared to the inverse of the 

.„), one 



minimal gap between adjacent frequencies {uiq 
has 



and 



k 



(24) 



(25) 



Using the orthogonality of the eigenvectors {|Aq)}, 
(Afc |A„i) = 5k,rn, one finally obtains the eigenvalue equa- 
tion of Cn : 



Cp|A,) 



(26) 



At equilibrium, the initial condition is thermalized and 
the energy is equally distributed among the modes. Each 
mode of frequency LOq and amplitude Uq carries an energy 
ma^^uj^ / 2 ~ kBT/2. Hence the eigenvalues Xq of Cp and 
the vibrational frequencies ujq of the dynamical matrix 
are related through 



A„ 



(27) 



A. Newtonian dynamics 

Starting from the Newton equations linearized around 
the reference state |r°), 

\Sr) +'D\Sr) =0, (19) 

where D = K/m is called the dynamical matrix, the 
eigenmodes of D, also called the vibrational modes of the 
system are defined by D|Aq) = u!q\Xq), where ujq are the 
vibrational frequencies. The solution of Ea. p9l) is given 
by: 

\Sr) ^e-'^'\Sr{0)), (20) 



B. Overdamped Langevin dynamics 

Starting from the Langevin equation for an over- 
damped linearized dynamics around the reference state 
|rO): 

\6rit))^~^\5r) + ^m), (28) 

where fi is the viscous damping and \ri{t)) is a white 
noise of amplitude F: {7]{t')\T]{t")) = T6{t' - t"), one 
introduces the operator C ~ whose eigenvalue 

equation is : C\Xq) = mujq / ^j,\Xq) . Tq ~ fi/nq = ^/(muj^) 
is the relaxation time of the system along the eigenmode 
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The solution of the Eq. (f28| can be written as 

1 rt 

\Sr) = e-^'\6r{0)) + - / e-^'-'-* '^T]{t'))dt' (29) 
M Jq 

Replacing this solution [25] in Eq. ([T]) one obtains an ex- 
pression with four terms, only two of which are not zero 
because the noise \ri{t)) is assumed not to be correlated 
with the position field \6r{0)). One ends up with: 



^MO)) = E,«JA,) and \rj{t)) = E.PMX,), with 
ctq = (A<j|i5r(0)) and Pq — {Xq\ri{t)), the first term of the 
previous equation, which we note Ti turns into: 



\ q,k 



\Xk){Xk 



e-^*|<5r(0)>('5r(0)|e-^*) + 



^^/g-£(t-t')|^(i'))(,;(t")|e-^(*-*": 



(30) where the last equality is justified by the fact that the 
components on each mode of the initial conditions are 
uncorrelated ( aqal) — aqa^^Sk.q- 

We now turn our attention to the second term in the 
Expressing |(5r(0)) and \rj{t)) in the eigenbasis {|Aq)} of Ea. (|3T|) . referred to as T2: 

I 



Mi- 



T2 



\{fdt" rdi'Ve-i/^'(*"*')/3,(^')|A,)(Afe|/3fe(t")e-i/-^(*-*"A 

yo Jo ^^fc / 



,-(t-t')/r,-(t-t")/rfc 



e - - ^ Pqit')Pkit"))\Xq){Xk\ 



r 



Now assuming that the components of the noise on the for particles in a newtonian fluid bath, but not necessar- 
modes are also uncorrelated ( Pq{t')f3k{t")) ~ r6q^kS{t' — ily for a large scale mechanical excitation. The generic 
t") one obtains for T2: equations of motion will then be: 



^2 = f dt'e-'('-'">/^^\X,){X,\ 

= 5Ey(i-^""'")l^^)(^^l 
k ^ 



Finally applying Cp on the eigenvector | Ag) one obtains 
the eigenvalue equation for Cp : 



Cp\Xq) — 



2 _ -2t/r, , £^ 

9 2/i2 ) ^ 2/i2 



(31) 



m\5f{t)) + ^L\5f{t)) = ~Vi\5r{t)) + \i^{t)) (32) 

We now expand in the modes of K , using the same 
notation as for the Langevin case, K|Aq) = Kq\Xq), and 
\5r{t)) ~ Yl,q'^q{^)W) ■ We define the colored noise in 
the basis of the modes as 

with {Vq{t)l^k{t'))=Tq5qk5{t~t'). (33) 

In the eigenbasis of the modes, the equations of motion 



C. Out of Equilibrium Stochastic Forcing 

To generalize on the two equilibrium cases of Newto- 
nian dynamics and Langevin dynamics, we now consider 
a non-equilibrium system with both inertia and damp- 
ing and a colored noise spectrum. We choose to stick 
to a Langevin type description of the dynamics, where 
the damping is linear and single particle, as appropriate 



,{t) + ^aq{t) + ^aq{t)=T^q{t). 



(34) 



We first solve the homogeneous equation, in the ab- 
sence of noise. The solutions are either oscillatory or 
overdamped depending on the relative importance of in- 
ertia and damping. 

For (/i/m)^ < AKq/m, i.e. the low damping limit, we 



19 



find two oscillating solutions 



For the oscillatory case, we obtain 



aq{t) = aq(0)e 2™ sin 



cos 



(35) 



corresponds to ncwtonian 



Clearly, the limit 11 /m 
dynamics. 

For (fj,/m)^ > AKq/m, i.e. the strong damping limit, 
we obtain two decaying solutions 



= = a^(0)cxp. 











m 





4^ 

m 



(36) 



It can be shown that in the limit m//^ — >■ 0, the '+'- 
solution corresponds to the Langevin homogeneous solu- 
tion aq{t) = Q;+(0)e~''9*/'^, while the ' — '-solution decays 
infinitely fast. 

To solve the inhomogeneous equation, we find a par- 
ticular solution a^{t) through the method of variation of 
constants based on the homogeneous solutions. In the 
oscillatory regime, we find 



<(t) 



/*!?iEl,-(t-O2sinf^(t-0 
Jo mU,„ \ 2 



(37) 



while for the damped regime, the particular solution is 
given by: 



'd<'!^e^(*-*')2sinh 
milq 



{t-t')] (38) 



Both solutions are similar in spirit to the particular so- 
lution of the Langevin equation: the source of continued 
motion in the system is the colored noise, and in the long 
time limit, only the motion due to the noise remains. 
From this point of view, the memory a non-dissipative 
system such as the Newtonian case retains of the initial 
conditions is singular. 

We can now calculate Cp by performing an ensemble 
average over initial conditions and the noise. In particu- 
lar, we assume that we can replace the noise correlations 
by their expectation value {'riq{t)Vk{t')) = ^qSqk5{t — t'). 
We also assume that noise and initial conditions do not 
cross-correlate, and that in ensemble-average, the initial 
conditions of different modes are independent of each 
other. Then only the diagonal terms in the modes re- 
main and, schematically, Cp is given by: 

Cp{t) = J2{aim + Y^iaim + ^«(o') (39) 



-fj,t/7n 



{{amr)sm^{nqt/2) 



{{a'^q{Q)f)cos\nqtl2) 
^iVqer'^ 

m 2lJLKq^l 



Tq 2r,e-^*/'" 



2lJLKq 



Qq sin ( flqt J cos ( riqt 



|A,)(A,| 
(40) 



For the damped case, we find 
Cp=J2 e-^*/" [(K(0))')e^^'' + ((a-(0))^)e- 



1 + e 



^ (e^^"* + e-^^'* - 2) 



TO A^KqQ^ ^ ^ ^ 

|A«)(A,| (41) 



which in the over-damped limit m/ fi — !■ 0, defaults to the 
result for Langevin dynamics if we assume equilibrium 
noise F^ = F, Vq. 

In both cases, in the long time limit we recover the an 
expression similar to that of the equilibrium result 



Cp\Xq) 



Z^Kq 



(42) 



where the initial conditions cease to matter. 



D. Sampling, Resolution and the Marcenko-Pastur 
theorem 



We derive relations which characterize the deviation 
induced by finite sampling and finite resolution on the 
spectral properties of Cp . For a system of N particles 
the positions of which are sampled independently n times 
with a resolution e, we will establish an adapted form 
of the Marcenko-Pastur theorem, which then allows to 
derive relations between the moments of the experimen- 
tally measured eigenvalue distribution and the moments 
of the eigenvalue distribution of the true correlation ma- 
trix. The sampling number r = Nd/n, is defined as the 
total number of degrees of freedom {NdY divided by the 
number of measurements n x Nd. 

In its most general form, the Marcenko-Pastur theo- 
rem [3^, [4^, [43, [mI is valid on the ensemble of Wishart 
random matrices, that is matrices constructed as Cr = 
X* X, where X is a n x p matrix which can be written 

1 /2 

'AS X = YCq where the elements of the n x p matrix 
Y are identically independently distributed (i.i.d), with 
mean 0, variance 1 and a finite fourth moment and Co is 
a, pxp positive definite matrix. In the following, we shall 
use the slightly more specialized derivation for Gaussian 
distributions found in Burda et al. [lH. We thus assume 
that Cr is the average of n samplings of XiXj, where the 
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variables {xi}, i = l...p are identically and independently 
distributed according to a multivariate normal distribu- 
tion: 



P{xi,...Xp) = [(27rfdet(Co)] 



-1/2 



exp 



1 



(43) 

In the limit n — >■ oo. at fixed p, one then clearly has from 
equation ([43|) limr^o C'r = Cq. 

For a given correlation matrix C, define the (negative) 
moment generating function of the eigenvalue spectrum 
P(A) 



This integral can be calculated by completing the square 
in the eigenbasis of Cq, and the resulting distribution is 
simply Gaussian with variance matrix Cq = Co + e^/. 

Then the correlation matrix which determines the 
probability distribution in our system is Co = Cq + e^/ 
and what is measured is Cr, such that lim„_i.oo Cr = Co. 
The Marcenko-Pastur theorem is still valid and we have 
Thr{zr) = 77io(2o)- We Can then relate the modified mo- 
ment generating function ?7io(2;o) to the real moment gen- 
erating function mo(zo) through a change of variables 
zq ~ zq ~ ■ The Stieltjes transform of the new eigen- 
value distribution becomes: 



k=l 



(44) 



where = (A''). Also define the Stieltjes transform of 
the eigenvalue distribution P(A) written as a sum over 
the poles of (C — zl)~-^ at the eigenvalues 

s(z) = -Tr((C-z/)-i) , zeC+. (45) 
p 



Both are related by 

m(z) 



-zs{z) — 1. 



(46) 



To sec this, note that in the eigenbasis of C, we can write 
the Laurent series of the right hand side as 



k=l ^ 1=1 



(47) 



Let us denote -Po(A) , mQ{z), so{z) the eigenvalue spec- 
trum, the generating function and the Stieljes transform 
of the true correlation matrix Co, and Pr{X), 771^(2;), 
Sr{z) the corresponding quantities for the experimental 
Cr- Then the Marcenko Pastur theorem stated in Burda 
et al. [51[ provides an explicit conformal transformation 
relating mr{z) to uiqIz): 



mr{zr) = 7710(2:0) where zq = 



1 + rmr{zr) 



(48) 



Burda et al. then derive relations between the moments 
ml. and through writing the Laurent series in Zr on 
both sides, and then equating the coefficients of the pow- 
ers oi z~^. 

We can adapt this approach to include an experimental 
gaussian noise term as follows. If Co is the true correla- 
tion matrix, let y.i ~ Sxi + be the particle fluctuations 
with the resolution-induced noise included, so that we 
have PiC^) = [27re2]-i/2exp(-^7e2). The probability 
distribution of a sum of two random variables is obtained 
through a convolution of their respective probability dis- 
tributions, so that we have 



P r 

P{y^,.■.,yp) = l[ J d^,P{iOP{yi-ii,-,yp-Cp)- 



(49) 



S{z) = iTr(C - zl)-^ = iTr(C - zl)-^ = s{z). (50) 

We can then derive the moment generating function: 

m{z) = - (zS(z)) - 1 = -(5s(z)) - 1 - e^siS). (51) 

The first two terms are nothing but 7*7(2), while one can 

show through a Laurent expansion that the last term is 

2 

given by ^[1 -1-771(2)]. The relation between the moment 
generating functions becomes then 



771(2) = 777(2) 



(52) 



from which we derive an adapted form of the Marcenko- 
Pastur theorem: 



777^(2^) = 7770(20) 
where zq 
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20' 



(53) 



1 -|- rmr{zr) 



Finally, this allows us to derive relations between the 
moments 777^? of the experimentally measured eigenvalue 
distribution and the moments of the eigenvalue distribu- 
tion of the shadow system uiq . This can be done through 
a Laurent expansion in Zr on both sides, together with a 
Taylor expansion in e^. After a considerable amount of 
algebra we find to second order in e^: 



771 J; 



777? 



777^ + 2777^6^ + r{?nl,f + 



(54) 



The central moments, i.e. the mean /i and the variance 
CT^, are more convenient and we find 



Mo 



(55) 



To this order, the second moment does not depend on 
the noise, however it comes in at higher orders. For a 
system with pure noise equation [SH] gives /7„ = and 
cr^ ~ re^ (note that for 7' — >■ 0, the second moment of 
the noise eigenvalue distribution vanishes as it should). 
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Since for large enough e, this has to be consistent with 
([55)) . we obtain: 

/ir = Mo + 

bl^al^ rul + re^ (56) 

To convert these relations into useful relations for the 
density of states one needs to perform the change of 
variable to oc A^^^^. Unfortunately, we could not find 
a straightforward and reliable way to calculate the mo- 
ments of this distribution and we had to rely mostly on 
numerical results. However in the limit r — > 0, the rela- 
tion C7o = Co -|- e^I directly translates to the eigenvalues, 
i.e. A^ = A^ + e^, and we recover equation [551 for r — Q. 



Here we can perform the change of variables explicitly, 
and we find to order (written in the most convenient 
mixture of direct and central moments): 

{alY = {alY-e^{u:%-{.^Uu.),) (57) 

It can be shown that the factor multiplying in the 
equation for is strictly positive, so that in the limit 
r — >■ 0, both mean and variance are reduced proportional 
to e^. 
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